Prony's method and the connected— moments expansion 



Francisco M. Fernandez 
INIFTA (UNLP, CCT La Plata-CONICET), Blvd. 113 y 64 S/N, 
Sucursal 4, Casilla de Correo 16, 1900 La Plata, Argentina 

We show that Prony's method provides the full solution to the nonlinear equations of the 
connected-moments expansion (CMX). Knowledge of all the parameters in the CMX ansatz is useful 
for the analysis of the convergence properties of the approach. Prony's method is also suitable for 
the calculation of the correlation function for simple quantum-mechanical models. 
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; I. INTRODUCTION 

Horn and Weinstein[l!| and Horn et alj2[ proposed the i-expansion for the calculation of the ground-state energy of 
' , quantum-mechanical models. It is the Taylor expansion about i = of a monotonically decreasing function E{t) that 
leads to the ground-state energy when t — >■ oo provided that the chosen reference function exhibits a nonzero overlap 
with the ground state. The coefficients of such Maclaurin series are known as connected moments or cummulants. 
^ Since the extrapolation of the t-expansion towards i — > cxi by means of Fade approximants[l[ Q did not appear 

^ ^ I ' to produce encouraging results, Cioslowski^S] proposed an exponential series and Stubbins 4] compared it with other 
extrapolation approaches. On matching the i-expansion and the exponential series at origin one has to solve a 



^ ' system of nonlinear equations. In order to bypass this problem Cioslowski^ developed a systematic algorithm for 
. . . 

obtaining just the parameter related to the energy by removing all the other variables in the nonlinear equations. The 
I resulting approach is known as connected moments expansion (CMX). From the properties of the Fade approximants 
' , KnowlesQ derived a compact an ellegant explicit expression for the approximants to the energy in terms of the 
^ I . connected moments. 

' Since the CMX approximants may exhibit singularities Q-I?! other authors proposed im pro ved approaches like the 
' alternate moments expansion (AMX)0, Q] and the generalized moments expansion(GMX) jlo|. Although these meth- 
ods may avoid the singularities in the CMX approximants they do not seem to improve the convergence properties of 
_ the sequences of approximants and commonly their results are not more accurate than the CMX ones. 

A recent analysis of the convergence properties of the CMX required the calculation of all the variables in the CMX 
ansatz For that reason, Amore and Fernandez pro posed a systematic method for the full solution of the CMX 
nonlinear equations motivated by an earlier discussionjisj of the connected-moments polynomial approach}l4j. 

Nonlinear equations like the CMX ones are an old problem in applied mathematics and were first solved by Prony[l^ 
and later Weiss and McDonough[l^ proposed an alternative approach based on Fade approximants. In fact, Frony's 
method is well known in numerical analysis [l^ and has been widely applied to a variety of problems in chemistry and 
engineering [isi liojl . Further inspection of the procedure developed by Amore and Fernandez 



12l | has revealed that it 



is closely related to Frony's method, although the main result derived by the former authors does not appear in the 



discussions of the latter method 
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The main purpose of this paper is to show the connection between Frony's method and the procedure developed 



2 



by Amore and Fernandez. In Sec. |TT]we show how to solve a particular set of nonlinear equations by means of Prony's 
method and derive the main result of Amore and Fernandezy . In SecM^e discuss the application of that method 
to the extrapolation of the generating functions for the moments and connected moments. In Sec. IIVI we test the 
general analytical results on simple models and show that they are useful for a discussion of the convergence properties 
of the CMX as well as for the approximate calculation of the autocorrelation function. 



II. PRONY'S METHOD 



In what follows we discuss the solution of the system of 2N equations 

N 



= fc=l,2,...,2iV 



(1) 



for the 2A'^ unknowns An and 6„, where Fk are known real numbers and s an integer. 
Prony's method is based on the construction of the polynomial [igI \v\ 

N N 

pip) = n - = Y.pi^^ = 1 

n=l j=0 

where the polynomial roots 6„, and thereby the polynomial coefficients Pj, are unknown. It follows from 

JV N N 

E p^+^p^ = E Y.p^bi, ^ 0, 1,2,..., N 

j=0 n=l j=0 

that the coefficients Pj are solutions to the linear system of equations 

E ^^+^p^ 

3=0 



(2) 



(3) 



^F^+N, I - 1,2,..., TV 



(4) 



If the matrix F with elements -Fi+j, i = 1, 2, . . . , A^, j = 0, 1, . . . , A^— 1 is nonsingular then the solution is unique. Once 
we have the polynomial coefficients pj we obtain its roots 6„ and then solve the resulting system of linear equations 
P]) (for example for fc = 1, 2, . . . , A^) for the remaining unknowns A^. 



The technique just outlined is known since 1795 



15j and is commonly called Prony's method|l6l. [l7 |. Weiss and 



McDonough[16j proposed an alternative way of solving the nonlinear equations by means of a partial fraction expansion 
of Pade approximants. and in what follows we describe a different strategy for obtaining the roots &„ developed by 



Amore and Fernandez 



12| . The starting point of this procedure is the homogeneous system of A^ linear equations 



AT 



E(^^ 



i+j 



bF, 



i+j-i 



ic,;=0, i = l,2,...,A^ 



(5) 



with A^ unknowns c^. There will be nontrivial solutions (c^ ^ 0) only if the determinant vanishes: 



Under such conditions we define 



N 
1=1 



(6) 



(7) 
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so that equations ([5]) become 7^ = b^j-i that leads to jj = 6^70, j = 1,2, . . . , N. Therefore, it follows from 

N N N N 

that the roots of the determinant ([5]) are exactly the roots of the polynomial p{h) of Prony's method. 

We have now at least three alternative procedures for solving the nonlinear equations ([1]). First, the original Prony's 
method that consists of solving the system of linear equations (|4]) for the coefficients of the polynomial p{h), then 
calculating its roots 6„, and finally solving the resulting system of linear equations ([1]) for the remainin g un knowns 



An- Second, the partial fraction expansion of the Pade approximants proposed by Weiss and McDonough[16|. Third, 
our approach that consists of obtaining the polynomial p{h) from the determinant ([S]). The choice of either of them 
is probably a matter of taste. We find our technique quite straightforward and it has revealed a most interesting 
connection between the exponential expansion of the generating function for the moments and the Rayleigh-Ritz 
method in the Krylov space fl^. 



3 



In what 



Prony's method was developed for fitting a series of exponential functions to given numerical data 
follows we apply it to a completely different problem so that the name Prony's method refers only to the strategy for 
solving the nonlinear equations ((T)). 



III. GENERATING FUNCTIONS FOR THE MOMENTS AND CONNECTED MOMENTS 

The general equations developed in the preceding section prove useful for the analysis of the CMX[3, 5]. We first 
consider the moments-generating function|]| 

Z(i) = (</)|e-*^|^) (9) 

where H is the Hamiltonian operator of the system and is an arbitrary reference state. The coefficients of its 
Maclaurin series 

^ i-ty 



^w-ES^/^.- (10) 



.=0 ' 

are the moments — |(/)). 

For simplicity we assume that the spectrum of H is discrete 

H \^,) = E, 1^,) , J = 0, 1, . . . (11) 

where Eq < Ei < E2 < ■ ■ ■, and without loss of generality we choose the eigenvectors of H to be orthonormal: 
ipj) = Sij. Under such conditions the generating function exhibits the exponential behaviour 

00 

Z(<) = ^|(0| V,)|'e-*^^- (12) 
j=o 

where the expansion coefficients |((/)| ipj)\^a.nd the energies Ej are unknown. We can calculate them approximately 
by means of the method developed in the preceding section and the ansatz 

JV-l 



where Aj and Wj are the approximations to ■0j)|^and Ej, respectively. If we require that the Maclaurin expansion 
of this ansatz yields the first moments fij exactly we are left with the nonlinear equations ([T]) with Fk = /ifc-i, s = — 1 
and Wj — 

It is has been proved that matching those expansions at origin is equivalent to the application of the Rayleigh-Ritz 
variational method in the Krylov space (RRK)}]^. Note that the determinant ([6]) becomes the secular determinant 
of the RRK. The RRK solutions 

JV-l 
i=0 

satisfy 

{<Pk\H\Vj)=W,{<Pk\^j) (15) 

and the application of Prony's method in the way just outlined yields Aj = \{(f>\ (pj)\^- 
The CMX is based on the generating function 



E(t) ^ (16) 



that is monotonically decreasing jll] and the coefficients of its Maclaurin series are the connected moments Ij: 



The recurrence relation 



°° (-fV 



h = A*! 



ij+i = Mj+i - E Ii+l^i.J-i, j = 1, 2, . . . (18) 

yields the connected moments (or cummulants) /; in terms of the moments u,-. 

n 

The main interest in E{t) is that it provides a size consistent approach to the ground-state energy[]| 

lim E{t) = Eq (19) 
when (01 V'o) 7^ 0. In order to carry out this extrapolation Cioslowskijsl proposed the exponential-series ansatz 

N 

E'^^\t)^AQ + Y,Aje'^^' (20) 

where the unknown parameters hj are supposed to be real and positive and Aq is the approximation to Eq. Matching 
this expression with the t-expansion (jl7p leads to the set of equations 

N 

Ao=h-Y,An (21) 

n=l 

and 

N 

Ik+i^Y.^^^n, k^l,2,...,2N (22) 
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We can solve the set of 2N equations (l22l) by means of the method of the preceding section with Fk = Ik+i and s = 0. 
In this case the exponential parameters bj, j = 1^2, . . . , N are the roots of the pseudo-secular determinant 

=0 (23) 

Once we have them we solve TV of the resulting 2N linear equations for the coefficients Aj, j = 1,2, ... , N, and 
then we obtain from equation pT|) . Cioslowskiy| and Knowles ^ developed remarkable strategies for obtaining 
^0 without calculating the other parameters explicitly. We have just shown that the explicit calculation of all the 
parameters in the exponential-series ansatz (1201) is quite straightforward. 



IV. ILLUSTRATIVE EXAMPLES 

In order to test the equations developed in the preceding section and show their usefulness in the analysis of the 
CMX we first consider an exactly solvable model: the dimensionless harmonic oscillator 



HQ 



As a trial function we choose one of the examples discussed in earlier papers 



= (x| </.) = (^x^ - e-2.V5 (25) 

for which -02)1 > |(0| '>pj)\, j 7^ 2. (The normalization factor is irrelevant to present purposes) One advantage of 
this example is that it is not difficult to obtain the generating function exactly 

, 121^3 + 189199^2 + 8180919U + 6561 

Fit) = ; t: , u — e 26) 

^' (81-u)(121m2 + 20198u + 81) ' ^ ' 



Note that limt^oo E{t) = Eq = 1 because {4>\ -ipo) 7^ in agreement with the discussion in Sec. 

By means of the general results of the preceding section we can easily calculate the generating function E{t) 
approximately in two ways: as f/^^^(t) = —Z'j^{t)/ZN{t) from equation (|13p and directly from equation (|20|) . Fig[T] 
shows that U'^^^t) approaches the exact generating function E{t) for all < t < 00 as increases. On the other 
hand, the approximate expression E^^'){t) given by Eq. (j20p is an unsatisfactory approximation to the exact E{t) as 
shown in Fig. [2] (except for sufficiently small values of t). This situation does not appear to improve as N increases. 
It is found that limt^oc E^^\t) = —00 for = 2,3 because in both cases there is one negative root bj associated 
to a negative coefficient Aj. For example, Ai « -0.0170, A2 « 0.147, ^3 « 0.000617, bi « -3.87, 62 ~ 4.04, and 
63 « 9.29 for = 3. Note that the curves in both figures are directly comparable in the sense that the calculation of 
either ZN+i{t) and E'^'^^t) requires exactly the same number of moments iij {2N + 1). The behaviour of E'^^^t) for 
larger N may be different; for example, for iV = 4, 5 there are two negative roots and limj^oo E^-^^{t) = +00. 

The anomalous behaviour of i?*^^^(<) is due to the poles of E{t) in the complex t-plane. Note that Eq. exhibits 
three real poles in the u-plane and the closest to the origin is located at u w —0.0040. It is reasonable to assume that 
the CMX ansatz E^^\t) should approach the Taylor series about u — for E{t) when the approach is successful. 
However, this expansion does not converge for u — I which is necessary for matching the i-expansion and exponential 



series at t = 0. This problem was discussed by Amore et alfllj by means of two-level models and we see it here again 
for the harmonic oscillator. Those authors have shown that the CMX approximants for the harmonic oscillator with 
the trial function (l^5t converge towards the second-excited state E2 ~ 5. From the approximate expressions E^^\t) 
with N ~ 1,2,3 we obtain Aq « 4.932,5.015,5.002, respectively, in agreement with those results. This example 
shows that the CMX approximants may converge to a meaningful result even when the ansatz (t) (on which the 
approach is based) is an unacceptable global approximation to E{t). 

Knowles|5[ argued that if the connected moments Ik are positive for all k then the bi are real and positive. Later, 
Massano et alQ] and Mancini et alP| suggested that Knowles' argument may not be valid and that it is the Hadamard 
determinants constructed from the Ik that should be positive. The example just analysed supports the latter conclusion 
because the first connected moments are all positive (« 5.13, 0.665, 2.12, 11.2, 40.0, 216, 979) but one of the roots 
is negative as shown above. We realize that the general formulas derived in Sec. |lT] are most useful for this kind of 
analysis. In fact, the calculation of the roots bj by means of present formulas is almost as easy as the calculation of 
the Hadamard determinants and the former provide a much clearer indication of the suitability of the ansatz E^-^"^ (t) 
as an approach for E{t). 

By means of the example discussed above we do not want to convey the wrong impression that the CMX is suitable 
for the calculation of excited-state energies because it is quite difficult to choose an appropriate trial function for 
this purpose (except when we can exploit the symmetry of the problem to make the trial function orthogonal to the 
ground-state [20j). The main interest in that exactly solvable problem is to show that we may predict an anomalous 
behaviour of the CMX by means of the roots bj . If they are not real and positive it is reasonable to suspect that the 
trial function may not be suitable and that it is convenient to choose another one. For example, if (f>{x) — all 
the roots bj are positive and the resulting £'(^'(t) is a remarkably good approximation to E{t) for all < t < 00 even 
at the low orders N = 1,2. 

Finally, we explore the possibility of calculating of the correlation function 

C{t) = Z{^t) = (^1 e'^'" !</>) (27) 

where |0) — |0(O)) is the initial state normalized to unity and \(f>{t)) = e^'*^ I'Pi^)) is the state at time t. In particular 
we concentrate on the real function |C(t)p. 

As a first example we consider the harmonic oscillator (p4|) and the trial function 



2X1/4 



(28) 



The exact correlation function is 



\C{t)\' = (29) 
^41 - 9 cos (4i) 

Fig. [3] shows resuhs for N = 2,3,4,5. There is no doubt that |Zjv(ii)| converges towards |C(f)| as increases. 
The approximation is satisfactory within the first period but the errors accumulate as t increases. This fact is not 
surprising because the approximate expression is constructed from the Maclaurin series for Z(it). 
The second example is the anharmonic oscillator 

H^-^+.^ (30) 
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and the same trial function ([28)) . Fig. |4]shows that the approximate results for = 2, 3, 4, 5 follow the same trend as 
in the case of the harmonic oscillator. At first sight it may be surprising that the convergence rate for the anharmonic 
oscillator appears to be greater than for the harmonic one. The reason is that the chosen trial state exhibits a greater 
overlap with the ground state of the former. Note that the magnitude of the coefficient Aq of Z^it) tells us that 
1(0 I'f/'o)!^ is approximately 0.943 for the harmonic oscillator and 0.981 for the anharmonic one. 
As a two-dimensional example we consider the anharmonic oscillator 

d2 cP 



where A > 0, and the trial function 



H^-^-^+x'+y^ + Xx'y'' (31) 

/r,\ 1/2 

0(x,y)=(^-j e-^'-y' (32) 

Fig. [S] shows results for A ~ 0.5 that look quite similar to those for the one-dimensional models. 

One does not expect that it may be possible to construct C{t) from its Taylor expansion about i = in all the 
cases. The problems discussed above are simple examples of single-well oscillators. In the case of the double-well 
oscillator V{x) — (x^ — 1)^ and the trial state localized in one of the wells 4'{x) = (f)^^^ e~'^^^' the approximants 
to |C(t)|^ do not appear to converge although the results for Z{t) look reasonable. If, on the other hand, we locate 
the initial state on top of the barier (for example, Eq. (f28| ) then the results are as accurate as those for the single 
wells. 



V. CONCLUSIONS 

The main goal of this paper is to show that Prony's method provides the full solution of the CMX equations in 
a straightforward and simple way. We have outlined the connection between Prony's method and the procedure 



developed recently by Amore and Fernandez 12|. Although both approaches are equivalent and Prony's one is known 
since long ago it seems that the pseudo-secular determinants obtained by Amore and Fernandez are not available 
elsewhere. 

In spite of the fact that the harmonic oscillator with the trial function (|25|) was chosen as an illustrative example in 
two earlier papers HQ 

we think that present discussion based on the exact generating function (I26p is clearer and 
provides more information about the behaviour of the CMX ansatz E'-^^t). In particular, we have shown that this 
simple problem clearly illustrates that Knowles' conclusion about the sign of the connected momentsj5|| is not valid 
as argued before by Massano et al[6] and Mancini et al[3|. 

It is clear that the roots bi provide a clear indication of the success of the CMX and the method developed by 
Amore and Fernandez jl^l , as well as Prony's method, make their calculation straightforward. Both procedures may 
be applied to the calculation of the autocorrelation function for simple quantum-mechanical models. If the approach 
converges then one obtains a reasonable analytic expression for the autocorrelation function within the first period of 
oscillation. 
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FIG. 1: Exact E{t) and approximants U^^^t) for the harmonic oscillator 
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FIG. 2: Exact E{t) and approximants _E'^^(t) for the harmonic oscillator 




FIG. 3: Exact |C(t)]'^ for the harmonic oscillator (line) and approximants \ZM{it)f with N = 2 (squares, blue), N — 3 (filled 
squares, red). A'' = 4 (circles, green) and N = 5 (filled circle, black). 
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FIG. 4: Approximants |Zjv(it)|^ for the anharmonic oscillator pO|) with N — 2 (squares, blue), A'^ = 3 (filled squares, red), 
= 4 (circles, green) and N = 5 (filled circle, black). 
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FIG. 5: Approximants |Zjv(it)!'^ for the anharmonic oscillator (|3ip with N — 2 (squares, blue). A" = 3 (filled squares, red), 
N = A (circles, green) and N = 5 (filled circle, black). 



